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Abstract. We consider the optimal control of quantum systems interacting 
non-linearly with an electromagnetic field. We propose new monotonically 
convergent algorithms to solve the optimal equations. The monotonic behavior 
of the algorithm is ensured by a non-standard choice of the cost which is not 
quadratic in the field. These algorithms can be constructed for pure and mixed- 
state quantum systems. The efficiency of the method is shown numerically on 
molecular orientation with a nonlinearity of order 3 in the field. Discrctizing 
the amplitude and the phase of the Fourier transform of the optimal field, we 
show that the optimal solution can be well-approximated by pulses that could 
be implemented experimentally. 



1. Introduction 

The control of quantum dynamics induced by an intense laser field continue to 
be a challenge to both experiment and theory [H [2 [3]. In this context, Optimal 
Control Theory (OCT) is an efficient tool for designing laser pulses able to control 
quantum processes [lllSlEllIllHlElIIOliniinillSlIIl]. By construction, the optimal 
field is the field steering a dynamical system from the initial state to a desired tar- 
get state and minimizing a cost functional which generally penalizes the energy or 
the duration of the field. Different methods have been developed to solve the opti- 
mal equations I15| . Among others, monotonically convergent iterative schemes 
proposed by Tannor et al. [16] and Rabitz et al. [H [171 [HI have been applied with 
success to a variety of physical and chemical processes [4j [HI [2Ql [211 [22 • These 
algorithms have the particularity to guarantee the increase of the cost functional 
at each step of iteration. In this paper, we will consider the Rabitz formulation 
of iterative algorithms [IT]- First introduced to treat pure-state quantum systems, 
these schemes have been extended and applied to mixed-state quantum systems, 
dissipative ones [21 [8] and non-markovian dynamics |23j . A majority of works has 
considered a linear interaction between the quantum system and the electromag- 
netic field. This linear interaction corresponds, for molecular systems, to the first 
order dipolar approximation (permanent dipole moment). Due to the intensity of 
the field or to the particular structure of the problem, some systems need to go 
beyond this approximation [211 [251 [Ml |27] . A typical example is given by the con- 
trol of molecular orientation and alignment of a linear molecule by non-resonant 
laser pulses [351 HSl [Ml HZj- When averaging over the rapid oscillations of the 
field, one observes that the permanent dipole moment plays no role in the control 
of the dynamics. In this case, molecular alignment and orientation are obtained 
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via the polarizability and the hyperpolarizabihty terms of the interaction Hamil- 
tonian (see [30] for information on the controUabiHty of these systems). From a 
methodofogical point of view, the natural question arises of whether one can apply 
monotonically convergent algorithms to such systems interacting non-linearly with 
the field. 

The goal of this work is to answer this question by proposing new monotonic 
algorithms when an arbitrary nonlinearity is considered. A key ingredient to en- 
sure the monotonic convergence of the algorithms is to consider a non-standard 
cost functional which instead of penalizing the intensity of the field, i.e., the square 
of the electric field penalizes a higher exponent which depends on the order of the 
non-linearity. Note that a similar question has been treated in |12j . A family of 
algorithms different from those proposed in this paper has been developed. In algo- 
rithms of [l^ , the cost is quadratic in the field and the control is decomposed into n 
components for a nonlinearity of order n. Thus, for each iteration of the algorithm, 
2n numerical resolutions of the time-dependent Schrodinger equation are required: 
n for the wave function and n for the Lagrange multiplier. On the contrary, in this 
work, we use only one component for the control field but at the price of modifying 
the cost functional. We construct monotonically convergent algorithms for pure 
and mixed-state quantum systems but they can be generalized straightforwardly to 
dissipative dynamics. We test the efficiency of these algorithms on the orientation 
dynamics of a linear molecule with non-linearity of order 3 corresponding to the 
hyperpolarizabihty terms of the molecule [271 . We use as target states the states 
which maximize the orientation in a finite-dimensional restriction of the Hilbert 
space. Several works have pointed out the role of these target states which both 
optimize the field-free orientation and its duration [3Iu32.,33. 34^. Promising results 
have been obtained both for pure and mixed-state quantum systems corresponding 
to zero and non-zero temperatures. 

Finally, we also analyze the structure of the Fourier transform of optimal control 
pulses. Our aim is to show that the optimal solutions can be well-approximated by 
pulses that could be implemented experimentally [351|351|371|3S] . Such pulses, tai- 
lored by genetic algorithms, have been successfully applied for experimentally and 
theoretically controlling different molecular processes [36l |37l |38l |39l |40l [41] . In 
the frequency domain, they are characterized by the fact that both the amplitude 
and the phase of the Fourier transform (but only for a finite number of frequencies 
equally distributed over a given frequency interval) are optimized [391 1401 141j . This 
choice corresponds to the types of pulses that can be implemented by liquid crystal 
pulse shapers. As an alternative, we use in this paper the results of our monotonic 
optimization algorithms to construct such pulses. Note that we do not adopt fil- 
tering techniques in the iterative algorithm, which have been proposed by several 
works (see [TTj and references therein). The idea consists generally in applying a 
filter to the control field at each iteration in order to satisfy spectral constraints. 
This filtering has the drawback that it does not generally yield a monotonic con- 
vergence of the algorithm. Instead, we propose to use a simpler solution. Starting 
from the optimal solution obtained by the monotonic algorithm, we discretize the 
phase and the amplitude of its Fourier transform into 640 points or less (640 points 
correspond to the number of pixels usually used in pulse-shaping experiments). 
From this discretization, we then construct a piecewise constant Fourier transform 
and a new time-dependent electric field by an inverse Fourier transform [40 [ I41| . 
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We finally compare the optimal result and the one obtained with the discretized 
field. We show that the difference between the two results is negligible when the 
structure of the optimal field is sufficiently simple or equivalently when the number 
of pixels is sufficiently large. 

This paper is organized as follows. We first present the model system in Sec. 
12.11 We determine in Sec. 12.21 the polynomial equation that must be satisfied 
by the optimal field. We then outline in Sec. 12.31 the principle of monotonically 
convergent algorithms for nonlinear interaction both in pure and mixed-state cases. 
A special attention is paid to the different choices and to the flexibility of the 
method. Generalizing the proofs of Refs. [17] and [22], we show the monotonic 
behavior of the algorithms. Section[3]is devoted to the application of these strategies 
to molecular orientation. The results are presented at T = K (Sec. 13. 2p and 
T ^ OK (Sec. 13. 5|) for the standard case (i.e. with a linear interaction term) and at 
T — OK for the averaging case (Sec. 13. 3p . We also propose an algorithm well-suited 
to the simultaneous optimization of two laser fields. An example is given by the 
non-resonant control of molecular orientation by two-color laser pulses [27]. We 
finally examine in Sec. 13.41 the structure of the Fourier transform of the optimal 
fields. 

2. Optimal Control Theory 

The goal of this section is to propose monotonically convergent algorithms suited 
to quantum systems interacting nonlinearly with the control field. To simplify the 
discussion, we consider the case of pure-state quantum systems. Following Ref. [8] 
and the formalism of super-operator, the proof can be straightforwardly extended 
to mixed-state quantum systems (see for that purpose Sec. 13. 5|) . Optimal control 
theory is invoked in order to maximize the projection onto a target state, but 
it could be equivalently defined for maximizing the expectation value of a given 
observable. The proof is a generalization of the standard proof for linear interaction 
[17] and of the proof given in Ref. ^22j . 

2.1. The model system. We consider a quantum system interacting with an 
electromagnetic field whose dynamics is governed by the following time-dependent 
Schrodinger equation 

(1) ^^^\m)-Hit)\m), 

which is written in units such that h = 1. The Hamiltonian H{t) of the system is 
given by 

(2) H{t) = Ho- pLE{t) - aE{tf - j3E{tf ■■■ , 

where Hq is the field-free Hamiltonian. The other terms describe the interaction 
between the system and the laser field E{t). This interaction is written as a polyno- 
mial expansion in E{t) whose coefficients are the operators fi, a, P ■ ■ ■ . For a linear 
molecule interacting with a linearly polarized laser field, the different operators jl, a 
and (3 are associated to the permanent dipole moment /io, the polarizability compo- 
nents a|| and a± and the hyperpolarizability components and (3± of the molecule 
[27] . These different molecular constants will be used in numerical computations of 
Sec. 121 
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2.2. Critical point. Let \4>o) and |(/)/) be the initial and the target states of the 
control. We denote by i/ the duration of the control. We define the optimal control 
theory through the following cost functional: 

(3) J=\{cbfmf))\^- r \E{tf-dt, 

Jo 

where n is a positive integer. The even exponent of the integrand and the choice 
A > ensure the negativity of the second term of Eq. ([3]). n is taken equal to 1 for 
a linear interaction but we will see that, in order to obtain monotonic algorithms, 
larger values of n have to be considered when the system interacts non-linearly 
with the field. A is a penalty factor which weights the importance of the laser 
fluence. Following [42], we will replace in Sec. [3] this constant by X/s{t), where 
s{t) = sin {iTt/tf), which penalizes more strongly the amplitude of the pulse at 
the beginning and at the end of the control. This allows to obtain more realistic 
optimal solutions. 

We introduce the augmented cost functional J which is defined through the 
adjoint state |x(*)) as follows 

J ^K0/iV'(i/))p - J^' \E{tf-dt - Hmmt]. 

where 3 denotes the imaginary part. The optimal electric field is solution of the 
equation 

which is a polynomial equation in E{t): 

(6) 2nAi;(t)2"-i + 2^[{i}j{tf)\dpf){x{t)\(i + 2aE{t) + ?,j3E{tf\ij{t))] = 0. 
The second term of Eq. ([6]) can be modified by using the fact that 

(7) |Wi)lxW)=0. 
The equation for the optimal field finally reads 

(8) 2nAi?(t)2"-i + 23[(VW|xW)(xWlA + 2ai?(t) + ?,mtf\m)] = 0- 

Setting the variations of J with respect to \'4'(t)) and \x{t)) to ensures that \^{t)) 
and |x(i)) satisfy the Schrodinger equation ([T]). To summarize, an extremum of J 
satisfies the equations 

{^,0 H{t)mt))^^ 



■dt 

\m) = \M 



(9) 

for the state \tp{t)) and 
(10) 

for the adjoint state |x(i)); the control field E{t) being solution of Eq. ([5]). 



\x{tf)) = 10/) 



MONOTONICALLY CONVERGENT OPTIMAL CONTROL THEORY OF QUANTUM SYSTEMS UNDER A NONLINEAR INTERA 



2.3. Monotonically convergent algorithm. We describe different iterative al- 
gorithms to solve the optimal equations of Sec. 12.21 To simplify the presentation 
of computations, we consider nonlinearity of order 3 and a cost which is quartic in 
the field (n = 2). 

At step fc > 1 of the algorithm, the system is described by the quadruplet 
(IV'fc(O): \Xk-i{t)),Ekit),Ek-i{t)) where |V'fc(i)) is the state of the system, \xk-i{t)) 
the adjoint state, Ek(t) and Ek-i{t) the electric fields associated respectively to 
the forward propagation of \ipk(t)) and to the backward propagation of |xfc-i(0)- 
lipkit)) and \xk~i{t)) are solutions of the following time-dependent Schrodinger 
equations 

(11) i^^\^k{t)) = HiEk)\Mt)) 
and 

(12) i^^\xk-i{t)) ^ H{Ek-i)\Xk-i{t)), 

where H[E{t)) = Hq- p,E{t) - aE{tf - j3E{tf . For \ipk{t)), we impose the initial 
condition |V'fe(0)) = |</>o) and for \xk-i{t)) the final condition \xk-i{t})) = \4'f)- 
The iteration is initiated by a trial electric field Eo(t). At step of the algorithm, we 
propagate forward the state \ipo{t)) with the electric field Eo{t). The cost functional 
Jfc at step k is defined by 

(13) Jk = \{<f>f\Mtf))\^ ~ r 

Jo 

The algorithm determines the quadruplet (lipk+iit)) ,\xk{t)) , Ek+i{t), Ek{t)) at 
step fc -|- 1 from the one at step fc. This is done by requiring that the variation 
AJ = Jk+i — Jk of the cost J from step fc to step fc -I- 1 is positive and that the 
limits (if they exist) of the sequences {Ek)k£T'! and {Ek)k£N are solutions of Eq. ([5]). 

For that purpose, we introduce the functions Pk+i{t) = \{xk{t)\tpk+i{t))\'^ and 
Pk+i{t) = \{xk+i(t)\ipk+i{t))\'^ ■ Differentiating with respect to time these two 
functions leads to 

(^)f/c+i(0 = fJ'k,k+i{Ek - Ek+i) + ak,k+i{El - i?fe+i) + l3k,k+i{El - i?|+i) 
for Pk+i and to 

—Pk+i{t) = ^I.k+l,k+l{Ek+l - Ek+i) + afe+i,fe+i(£'fc+i - E'fc+i) 

(15) +[ik+i,k+i{Elj^i- EIj^{) 

for Pk+i- In Eqs. ([Till and (fT5|) . we have introduced the notation 

(16) AkM' - 2^[{^k'{t)\xk{t)){xk{t)\A\^k'm 

for a given observable A. The functions P^ and i\. fulfill by definition the following 
relations: 

Pk+i{tf) = mu^k+iitfW = Pk+i{tf) 

(17) Pfc+i(0) = |(xfc+i(O)l0o)|' = Pk+2m, 
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and a direct integration gives 

(18) p^(^tf)^P,iO) + J^'^dt. 

The variation A J is given by 

AJ 419^+1 -Jk^ l(0/IV'fc+i(i/))l' - \{<l^f\Mil))? - / >^[Ek+i{tf - Ek{tf]dt. 

Jo 

Using the fact that |(0/|V'fe+i(t/))P - \{(t>f\Mtf))? = Pk+iitf) - Pkitf) and Eqs. 
(flT)) . one deduces that A J — Pi + P2 where 

Pi^- r ^[Et+i - + 

Jo 

(20) / [{Ek-Ek+i)fikM+i + {El-El+,)ak,k+i + {El-El+,)(3k,k+i]dt 
Jo 

and 

P2^ r X[Et-Et]- 
Jo 

(21) / [{Ek - Ek)iik.k + {El - El)ak.k + (El - El)h.k\dt. 
Jo 

To ensure the monotonic behavior of the algorithm, we choose the fields Ek+i and 
Ek such that the integrals Pi and P2 are positive. A sufficient condition is to 
impose that the two integrands Vi and V2 associated to Pi and P2 are positive 
[22j . To be more precise, we first determine Ek from E^ such that P2 is positive 
and then we determine E^+i from Ek such that Pi is positive. \'ipk+i) and \xk) are 
computed from a forward and a backward propagation with the fields Ek+i and Ek- 

Starting from these conditions, we introduce two algorithms. 
Algorithm I: 

Vi and V2 are respectively viewed as functions of Ek+i and Ek- Ek+i and Ek 
are defined as the control fields which maximize Vi and V2- The maxima of these 
polynomials are positive since Vi{Ek) ~ and V2{Ek) — 0. As already mentioned, 
we first determine for each time t the maximum of 'P2 and then the one of Vi. 
The integer n of the cost is chosen sufficiently large to ensure that the fields which 
maximize Vi and 1^2 are finite. This means that we choose n such that the terms 
— A-B^"]^ and — Ai?^" are the monomials of higher degree in Vi and V2- We then 
have: 

lim Pi(Sfc+i) = 

-Efc+i^±oo 

^ lim V2{Ek) = 

»±oo 

which satisfies the requirement. For nonlinearity of order 3, the choice rt = 2 is 
sufficient. Finally, if we assume that the algorithm converges then we can check 
that the limits of the sequences {Ek)keN and {Ek)keN are solutions of Eq. ([8]). For 
that, we respectively differentiate Vi and 7^2 with respect to Ek+i and Ek- We 
next replace in the derivatives of Vi and V2, Ek+i, Ek and Ek by E, \ijjk+i), iV'fc) 
by \tjj) and \xk) by \x)- It is then straightforward to see that the limit E{t) satisfies 
the optimal equation ([8]). 
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Algorithm II: 

We first write Vi and V2 as follows: 

(22) Vi = (Eu+i - Ek)[-X{El+^ + El^^Ek + E^+iEl + E^ 
-/i/c,fe+i - ak,k+iiEk+i + Ek) - Pk,k+i{El^i + EkEk+i + Ef)], 

and 

(23) V2 = [Ek - Ek)[-\{El + ElEk + EI^ + El) 
-Mfc^fe - ak,k{Ek - Ek) - Pk,k{El + EkEk + El)]. 

We then introduce two positive constants rji and 772 by setting: 

Ek+i ~ Ek ^ r]i[~\{El_^^ + El^^Ek + Ek+iEl + El) 

(24) — — (Xk^k+iiEk+i + Ek) — Pk,k+iiEl^i + EkEk+i + Ef.)], 
and 

Ek-Ek^ m[~KEl + ElEk + ElEk + El) 

(25) -^lk^k ~ ak,k{Ek + Ek) - h,k{El + EkEk + El)]. 

Equations. ([24|) and ((25|l are viewed respectively as equations in Ek+i and ii^/j. Ek+i 
and are defined as one of tlie solutions of these two equations. By definition of 
the constants ryi and 772, the values of Vi and V2 for these fields are positive. The 
integer n is chosen sufficiently large to ensure that Eqs. ((24| and ((25l) always have a 
real solution respectively in Ek+i and Ek. For nonlinearity of order 3, it is sufficient 
to take n = 2. When Eqs. ([24|) and ([25|) have more than one real solution at time 
t, we numerically choose the solution that is closest to the one at time t — dt (for 
the forward propagation) ot t + dt (for the backward propagation) . The processus 
is initiated by imposing that Ek{0) — and Ek{tf) = 0. This allows one to obtain 
smooth optimal fields without discontinuity. As for the algorithm I, we can check 
that the limits of the sequences {Ek)keN and {Ek)k£N satisfy Eq. ([5]). This can be 
done by replacing Ek and Ek+i by E in Eqs. ((24l) and ((25|l . 

In the two cases, the structure of the algorithms can be summarized as follows. 
At step fc + 1, we propagate backward in time the adjoint state \xk) with the field 
Ek determined from P2. We then compute the forward evolution of l^'k+i) from 
|0o)- For this second propagation, we use the field Ek+i defined from Pi. Note that 
a simpler solution which gives a slower convergence consists in choosing Ek — Ek, 
i.e., to propagate \ipk) and \xk) with the same field. 

3. Control of molecular orientation 

3.1. Introduction. In this section, we investigate the control of orientation dy- 
namics of a diatomic molecule driven by an electromagnetic field [28l |29]. This 
control is taken as a prototype to test the efficiency of the algorithm. The applica- 
tion of OCT to molecular alignment and orientation is relatively recent p2t \12\ [43] . 
One of the main results of Ref. [52] is that the optimal oriented state (see below for 
a definition) is reached by rotational ladder climbing, i.e., by successive rotational 
excitations. The corresponding optimal pulse is however very long, of the order of 
20 rotational periods, which could be problematic for practical applications. We 
consider shorter durations in this paper of the order of the rotational period Tper- 
We have chosen t/ = Tpe,. but other durations can be considered. Note that for 
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controls much shorter than Tper, the optimal solution is very close to the kick mech- 
anism largely explored using the sudden- impact model [44]. The CO molecule is 
taken as an example. The units used are atomic units unless otherwise specified. 

The molecule is described in a rigid-rotor approximation interacting with a lin- 
early polarized laser pulse nonresonant with vibronic frequencies. In this case, the 
Hamiltonian H can be written as follows [45l [27] 



where B and /xq are the rotational constant and the permanent dipole moment. 
a||, a±, /?|| and P± are respectively the polarizability and the hyperpolarizability 
components of the molecule. The labels || and ± indicate the components parallel 
and perpendicular to the internuclear axis. For the CO molecule, we have chosen 
the following numerical values B = 1.9313 cm"-'^ and /Zq = 0.044, a|| — 15.65, 
a±_ = 11.73, /3|| = 28.35 and l3±_ = 6.64 in atomic units [H[17]. is the angular 
momentum operator and 9 the angle between the direction of the molecular axis 
and the polarization vector. A basis of the Hilbert space is given by the spherical 
harmonics \ j, m) with j > and — j < m < j. 

3.2. Zero rotational temperature. In this section, we consider the limit of zero 
rotational temperature. We recall that the expectation value (cos 9) is usually taken 
as a quantitative measure of orientation [28l [29]. Here, we replace this measure 
by the projection onto a target state \(j)f). We consider target states recently 
introduced for the orientation which both maximize the field-free orientation and 
its duration [32l[3T]. To construct this target state, we restrict the Hilbert space 
to a finite-dimensional one defined by a maximum value of j denoted jopt- For 
CO, we have chosen jopt = 4 which leads to a maximum of (cos 9) of the order of 
0.9. In this reduced Hilbert space, the operator cos9 has a non-degenerate discrete 
spectrum. The target state \(j)f) is then defined as the eigenvector of cos of highest 
eigenvalue. The initial state is the state jO, 0). We also recall that the projection m 
of the angular momentum j on the field polarization axis is a conserved quantum 
number due to cylindrical symmetry. 

We now apply the monotonically convergent algorithms I and II. The results of 
the computations are presented in Figs. [TJ [21 [21 and [Jj We have used the simplified 
algorithms by assuming that Ek ~ Ek ■ Figures [TJ [3| and [4| correspond respectively 
to the algorithm II for n = I, the algorithm II for n — 2 and the algorithm I for 
n ~ 2. For the algorithm II, we can choose n — 1 since we have checked that for 
this value Eqs. p4)l and ([25]) always have a real solution. Note that this latter 
observation depends on the values of A and 77 considered. Numerical values are 
taken to be A = 0.05, 77 = 1 for Fig. [T] (A corresponds to the maximum value of 
A(t)), A = 6.05 X 10^,77 = 1 for Fig. [Hjand A = 12 x 10^ for Fig. [D The difference 
in the values of A is due to the form of the cost which is either quadratic or quartic 
in the field. We have checked that the value of rj is not relevant even if the value 
of A has to be adjusted with respect to the one of 77. The trial fields are displayed 
in Fig. [1^, [3^ and [4^. The trial field is a gaussian pulse of intensity of the order 
of 1 TW/cm^. In order to obtain realistic electric fields, the value of A has been 
chosen so that the energy of the optimal field be lower than two times the energy 
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(b) 




Figure 1. Plot as a function of the adimensional time t/Tper of 
(a) the optimal field (solid line) and the initial trial field (dashed 
line), (b) the expectation value (cos 6*) and (c) the projection onto 
the target state 10/). The abbreviation a.u. corresponds to atomic 
units. 




(a) 




(b) 



10 20 30 40 50 

number of iterations 



Figure 2. Plot as a function of the number of iterations of (a) 
the adimensional cost J defined by Eq. for n = 1 (a cost 

quadratic in the field) and (b) the projection onto the target state 
at time tf. 



of the gaussian pulse. Note also that for E ~ 5.10^^a.u. (which corresponds to the 
typical amplitude of the optimal field), we have /xq — a_\_E which shows that the 
polarizability terms are not negligible in the dynamics. In each case, very good 
results are obtained with a final projection |(V'(^/)|^^'/)P larger than 0.99 except 
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Figure 4. Same as Fig. [T]but for the algorithm I and n — 2. 
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for the algorithm I where | (V'(i/)I0/) P is of the order of 0.98. Figure [2] illustrates 
the convergence properties of the algorithm II which are satisfactory since after 
30 iterations, we obtain a projection close to 0.98. A similar behavior has been 
observed in the other cases. A comparison of Figs. [T^, [3^ and|4^ shows that the 
optimal field for n = 2 has sharper variations than for n = 1 for both algorithms. We 
have also observed that these sharper variations can induce numerical instabilities 
and high frequency oscillations in the optimal field. This point is discussed in Sec. 
13.41 where we show how to remove the parasite oscillations with a band-pass filter. 
In practice, it has been found that small values of the exponent n generally produce 
smoother optimal fields. 

3.3. Non-resonant two-color laser fields. We continue to consider a zero ro- 
tational temperature but we assume now that the molecule interacts with a non- 
resonant two-color laser field [24l [27] of the form 



(27) E{t) = Ei{t) cos(wi) + E2{t) cos{2ujt). 



H^BJ^ - -[(«!! - a_L) cos^ 9 + ai]{Ei{tf + E2{tf) 



After averaging over the rapid oscillations of the field, the Hamiltonian H of the 
system becomes 

(28) -3P±)cos^e + 3/3_LCOsd]Ei{tfE2{t). 

8 

The interest of this model is due to the absence of linear term in the interaction 
which enhances the difficulty of the control. Two cases can be considered according 
to the respective values of Ei and E2- If Ei = E2, we can use the standard 
algorithm presented in Sec. 12.31 whereas for Ei E2, the algorithm has to be 
slightly generalized. These two problems are respectively analyzed in Sees. 13.3.11 
and [3X1 

3.3.1. The case Ei ^ i?2. We first generalized the algorithm of 12.31 to the case of 
two control fields. We assume that the cost is quadratic in the field. 

We introduce the augmented cost J and we determine the critical points with 
respect to Ei and i?2. We have: 

J=\{<t>fmtf))\' - X[E,itf + E2{tf]dt 







d 



(29) -2^mfMf)J^ (^^(t)\{i--H)m))dt], 

and we compute the variational derivatives 6J/6E1 and SJ/SE2 which are equal to 
zero for a critical point. We then obtained the following system of equations 

XEi + 2&E1 + 2PE1E2 - 
^ ' \E2 + 2&E2 + (3El = 

which are satisfied by the optimal fields Ei and E2- We have used in Eqs. ((30| the 
notations 

a = 2S[(^(i)|xW)(xWl3((a|| " a±) cos^ 9 + a^)^))] 
^ i3^2^[(m\x{t))(x{t)\lm~3P^)cos^9 + 3(3^cos9]\m)]- 
We solve the optimal equations by a monotonically convergent algorithm. The 
proof of monotonicity follows closely the lines of proof in Sec. 12.31 We use the 
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same notations, with for instance Ei k the field Ei at iteration k. To simpUfy the 
computations, we take equal the fields and Ek for the forward and the backward 
propagations. We compute AJ — Jk+i — Jk- We obtain the following expressions 
for the polynomials Vi and 1^2'- 

(32) Vi = {Ei^k+i - Ei^k)[{Ei,k+i + Ei^k)ak^k+i + 
{Ei^k+i + Ei^k)E2,k+iPk,k+i — A(i?i.fc+i + Ei^k)] 

and 

(33) V2 — {E2^k+i — E2,k)[{E2,k+i + E2,k)o!k,k+i + 

EfkPk,k+l — A(£^2,fc+1 + E2,k)] 

where 

.m,fe+i = 9[(V'fc+i(i)IXfc(i))(Xfc(i)l?((a|| + a±)cos2 0|i/'fc+i(t))] 

W,fc+i = S5[(V'fe+i(i)IXfe(t))(Xfc(i)l|[(/3|t -3/5i)cos='0 + 3/3iCos0]|V^,+i(O)]. 

■Pi and V2 are respectively viewed as polynomials in £'i,fc+i and E2.k+i- We first use 
7^2 to determine the field i?2.fc+i by the algorithm I or II and then using this solution, 
we compute Ei^k+i from Vi- We also check that if the algorithm converges then 
the solutions given by the algorithm correspond to the extremal solutions defined 
by Eqs. ((30)) . This can be done by replacing Ei^k+i and Ei^k by Ei and E2^k+i and 
E2,k by £'2. The optimal fields are then zeros of the derivatives of Vi and V2 with 
respect to Ei^k+i and i?2,fc+i- Figure [5] displays the results we have obtained with 
the algorithm II for 71 = 2 and A = l.We have chosen two different trial fields in 
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Figure 6. Same as Fig. [T]but for the averaging case with Ei — £2- 

order to generate two different optimal fields Ei and £2- With the same trial field 
for the fields Ei and E2, the algorithm leads to two solutions which are very close 
to each other. Larger values of electric fields have been used due to the absence of 
linear interaction term in the Hamiltonian. 

A remarkable characteristic of the optimal fields is the fact that Ei vanishes 
for t > 0.2 X Tper- This means that the dissymmetry producing the orientation 
(dissymmetry due to the term in E\E2 in the Hamiltonian) only acts during this 
duration. This provides a non-intuitive and new method to produce orientation 
using a long laser field E2 and a short laser field Ei . 

3.3.2. The case Ei = i?2. We use the monotonic algorithm II proposed in Sec. 
m Figure [H] illustrates the different results. They have been obtained for n = 2 
and A = 5. The trial field is a gaussian pulse whose duration corresponds to the 
rotational period. From the equations of the algorithm, it is straightforward to see 
that the algorithm cannot generate an optimal field different from zero at time t 
if the trial field is zero at that time. This is simply due to the absence of linear 
interaction term in the Hamiltonian. The algorithm only modifies the envelope of 
the trial field whose choice is therefore crucial. 

3.4. Analysis of the Fourier spectrum. We analyze in this section the Fourier 
transforms of the optimal solutions. Our goal is to show that optimal solutions 
determined in Sec. I2.1l can be well-approximated by solutions that could be imple- 
mented experimentally. We consider experiments coupled with genetic algorithms 
optimizing the phase and the amplitude of the Fourier transform of a finite num- 
ber of frequency components. The discretization is done over a frequency interval 
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Figure 7. Plot as a function of time t/Tper of (a) the optimal 
(dashed line) and of the approximate fields (solid line) and (b) the 
projection onto the target state |^/). 




(a) 





Figure 8. Plot as a function of the frequency ly of the module 
square of the Fourier transform of the optimal field (a) and of its 
piecewise constant approximation (b). The Fourier transform has 
been discretized over the interval [—1^4, 



chosen with respect to the quantum transition frequencies involved in the control 
(see below). Note that we do not take into account, in this paper, technology con- 
straints for the choice of this frequency interval. Standard pulse shapers usually 
work with optical frequencies of the order of 800 nm. With such a technology, only 
non-resonant laser fields of Sec. 13.31 could be experimentally implemented. 

Following [40l |4T] , we assume that the solution obtained by genetic algorithms 
is a piecewise constant function in frequency both in amplitude and in phase. We 
have chosen 640 frequencies or less to discretize the optimal field. By an inverse 
Fourier transform, we then determine a new time-dependent electric field. Figure 
[7] presents the results obtained with the optimal pulse of Fig. [1] Using only 128 
frequencies, we show that the final projection obtained by the optimal pulse and its 
approximation are very close to each other. For 256 frequencies, the difference is 
neghgible and cannot be distinguished at the resolution of the plots. Figures [8] and 
[9] give informations on the Fourier transform of the optimal pulse. One introduces 
the rotational frequencies Vj+i given by: 

(35) i/j+i = E,+i - E, = 2B{j + 1) 
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frequency (v) 

Figure 9. Same as Fig. but as a function of the frequency v. 
Vertical lines indicate the positions of the frequencies i/i, 1^3 
and V4. 

where Ej is the energy of the state |j, 0). The target state being associated to 
jopt = 4, only five rotational states from j = to = 4 have to be populated by the 
control field. It is thus natural to discretize the Fourier transform over the interval 
[—1^4:, 1^4]. Higher frequencies do not contribute to reach the target state. A similar 
behavior has been observed for the other optimal solutions. Another example is 
given by Fig. 1101 The corresponding optimal solution obtained by the algorithm 
I presents rapid unwanted oscillations. To obtain a smooth solution displayed in 
Fig. llOb . we filter this optimal pulse in the frequency domain. The bandwidth 
of the filter is chosen to cut off frequencies higher than 1^4 which produce rapid 
oscillations. As in the first case, we also observe that the discretization does not 
significantly modify the final result. 

3.5. Finite rotational temperature. We investigate the temperature effects on 
the optimal solutions. The system is described by a density matrix p whose dynam- 
ics is governed by the von Neumann equation. The initial density operator p{0) is 
the equilibrium density operator at temperature T which can be written 

1 +00 j 

(36) piO) = ^ E E e-^-'(^+i)/(^--^) \j, m) (j, m| 

where ks is the Boltzmann constant and Z the partition function. The objective 
of the control is to maximize the projection of p(tf) onto a target state Popt- We 
consider here the target state introduced in Ref. [33] which is unitarily equivalent 
to the initial mixed state p{ti) and optimizes both the orientation and its duration. 
We refer the reader to Ref. [33] for the complete construction of Popt and for proofs 
of its attainability by unitary controls. Note that the definition of Popt depends 
on the polarization used. We consider here the optimum for a linear polarization. 
Popt can be defined as follows. The first step consists in reducing the dimension 
of the Hilbert space to a finite one Ti^J'pt) where jopt is the highest j for which 
the corresponding rotational levels are significantly populated. The dimension of 
this space depends on the temperature and on the intensity of the field used. For 
CO and T = 1, 5 and 10 K, we have chosen jopt = 4. We denote by H^""*' the 
subspace of Ti.^^""*^ associated to a given value of m. The target state pi!"/*' of the 
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Figure 10. Plot as a function of time t/Tper of (a) the optimal 
(solid line) and (b) the approximate fields (dashed line) and (c) 
the projection onto the target state \4>f). 



control, which therefore depends on the choice oi jopt, is given by 

m=jopt jopt-\ni\ + l 

(37) piirr'= E E '^rvr'xxn, 

where the Wj;™'''s are the eigenvalues of p{ti) restricted to Hm'""'^ and ordered. The 

vectors Ixi™^) ^^re the eigenvectors of the restriction of the operator cos 

The vectors Ixi™^) ai"© also ordered according to the values of the corresponding 
eigenvalues. 

We have used the algorithm II with n = 2 to determine the optimal solutions. 
We denote by xi^) the adjoint density matrix state. In superoperator notations, 
the structure of the algorithm is very similar to the one for pure states. The cost 
functional J is given by 

(38) J = \{{PopMtf)))f - r \E{tfdt 

Jo 

and the augmented cost functional J reads 

J= \{{PopMtf)))f-2m{p{tf)\p,pt)) f\{x{t)\i.^^^-H)\p{t)))dt] 







(39) - / ' \E{tfdt 

Jo 



MONOTONICALLY CONVERGENT OPTIMAL CONTROL THEORY OF QUANTUM SYSTEMS UNDER A NONLINEAR INTERA 




Figure 11. Plot as a function of the adimensional time t/Tp^r 
of (a) the optimal electric field and (b) the projection onto the 
target state Popt- Solid, dashed and dot-dashed lines correspond 
respectively to T = 1 K, T = 5 K and T = 10 K. 



where {{x\p)) = Tr[xV] and {{x\M\p)) = Tr[x^[M,p]] for a given observable M . 
p{t) and x(<), which satisfy the von Neumann equation, are propagated forward 
and backwards with initial condition p(0) — pQ and final condition x(i/) = Popt- 
Numerical parameters are respectively taken to be A = 6.10^, A = 65.10^ and 
A = 90 for T = 1 K, T = 5 K and T = 10 K. 77 is equal to 1 in aU the cases. The 
values of A are chosen so that the total energy of the field stays approximatively 
constant when the temperature is increased. The trial field is the same for the three 
cases considered. 

Figure [TT] presents the results obtained for three different temperatures. We 
observe that the structures of the optimal fields and of the projection onto the target 
state as a function of time are very different. As expected, we note a decrease of the 
final projection with increasing temperature. This computation allows us to show 
the robustness with respect to temperature of the optimal fields. For T = 10 K, 
we still obtain an efficient field since the final projection is of the order of 0.9. We 
also point out sharp variations of the optimal fields due to the use of a cost that is 
quartic in the control field. These variations do no affect the temporal evolution of 
the projection. 



4. Summary 

We have presented a new family of monotonically convergent algorithms for the 
computation of the optimal control of a quantum system interacting non-linearly 
with the laser field. One key for the convergence of these algorithms is to consider 
costs which are not quadratic in the field. In comparison with algorithms of Ref. 
[12j . this allows to consider only one wave function and one adjoint state per itera- 
tion of the algorithm whatever the nonlinearity used. This is thus less demanding 
from a numerical point of view especially when the degree of the nonlinearity is 
important. As a prospect, an open question in this field is the applicability of the 
present method to more complicated systems involving, for instance, non-markovian 
dynamics or a time-dependent target state |48j . Special attention has to be paid 
to the convergence properties and to the stability of the method, especially when 
the cost is quartic in the control field. 
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